Numerical simulations of stiff fluid gravitational singularities 
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Numerical simulations of the approach to the singularity in spacetimes with stiff fluid matter are 
presented here. The spacetimes examined have no symmetries and can be regarded as representing 
the general behavior of singularities in the presence of such matter. It is found that the singularity 
is spacelike and that as it is approached, the spacetime dynamics becomes local and non-oscillatory. 
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I. INTRODUCTION 

A longstanding problem in general relativity has been 
to find the general behavior of singularities. Several re- 
sults, both analytical[l| and numer ical@ have been ob- 
tained. Though most of the results are for the case where 
the spacetimes have one or more symmetries, recent work 
has been done on the general case where there are no 
symmetries. 0- EJ There is a longstanding conjecture^ 
due to Belinski, Lifschitz and Khalatnikov (BKL) that 
states that the generic singularity is spacelike and local. 
This conjecture has been reformulated and put more pre- 
cisely by Uggla et al. |7f The type of local dynamics con- 
jectured by BKL depends on the type of matter. For vac- 
uum and for many types of matter, the BKL conjecture 
is that the local dynamics is oscillatory, corresponding to 
the dynamics of a Bianchi type IX spacetime. However, 
for stiff fluid (i.e. fluid with pressure equal to energy 
density) the BKL conjecture is that the local dynamics is 
asymptotically velocity term dominated corresponding to 
the dynamics of a Bianchi type I spacetime. The vacuum 
version of the BKL conjecture has been supported by the 
numerical simulations of p| which show local and oscilla- 
tory dynamics in vacuum spacetimes with no symmetry. 
The stiff fluid version of the BKL conjecture has been 
supported by the theorem of Andersson and RendalLSj 
which shows the local existence in a neighborhood of the 
singularity of solutions of the Einstein equations with stiff 
fluid matter with the expected asymptotic behavior and 
with enough degrees of freedom to be the generic solu- 
tions. What is not known is whether generic stiff fluid 
initial data evolves to a solution of the Andersson and 
Rendall class. 

To address this issue, we perform numerical simula- 
tions of the approach to the singularity for stiff fluid 
matter with no symmetries. Our methods are those of Q 
using the system of . Section II presents the equations 
and numerical methods used. The results are given in 
section III and conclusions in section IV. 



II. EQUATIONS 

The system evolved here is essentially that of 
reference^ but specialized to the stiff fluid case and with 
a slightly different choice of gauge. Here the spacetime 
is described in terms of a coordinate system (t, x 1 ) and a 
tetrad (eo,e a ) where both the spatial coordinate index 
i and the spatial tetrad index a go from 1 to 3. It is 
assumed that eo is hypersurface orthogonal and that the 
relation between tetrad and coordinates is of the form 
e = N~ 1 d t and e a — e a l di Here N is the lapse and 
we have chosen the shift to be zero. We choose the spa- 
tial frame {e a } to be Fermi propagated along the integral 
curves of eo. The commutators of the tetrad components 
are decomposed as follows: 



[e ,e Q ] = u a e - (H5 a ^ + cr Q p )e /3 



(1) 
(2) 



where n a P is symmetric, and a a ^ is symmetric and trace 
free. Square brackets denote the antisymmetric part of a 
tensor. 

Define u a = eQ a and h ab = g ab + u a u b , that is u a is the 
timelike vector of the tetrad and h ab is the spatial metric 
corresponding to the choice of u a as the time direction. 
Then the stress-energy tensor can be decomposed as 



Tab = fJ.U a U b + 2q {a U b ) + ph ab + TT ab 



(3) 



where q a and ir ab are orthogonal to u a and where ir ao 
is symmetric and trace-free. Round brackets denote the 
symmetric part of a tensor. 

Scale invariant variables are defined as follows: 

{d ,d a } = {e ,e a }/H (4) 

{E a \Z a/3l A a ,N af3 } = {e a i ,a aP ,a a ,n aP }/H (5) 

q + l = -d liiH (6) 

r a = -d a \nH (7) 

{Q, P, Q a , n Q/3 } = q a , n af3 } /(3H 2 ) (8) 

The matter variables are not all independent, because 
we assume that the stress-energy is that of a stiff fluid 
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T ab = jl{2u a U b + g ab ) 



(9) 
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Here jl is the rest frame energy density of the fluid and 
u a is the fluid four-velocity, which can be decomposed as 
u a = T(u a + v a ) where v a is orthogonal to u a . Compari- 
son of equations © and © yields 



n 



Q - — v 

2n 



a/3 



G- 



-V <a Vj3> 



(10) 

(11) 
(12) 



Here v 2 — v a v a and G+ = 1 + v 2 and angle brackets 
denote the symmetric trace-free part of a tensor. Thus, 
all scale invariant matter variables can be expressed in 
terms of fl and v a . 

Finally choose the lapse to be N = H~ l . The relation 
between scale invariant frame derivatives and coordinate 
derivatives is 6q — dt and d a = E a l di. From the Einstein 
field equations and the conservation of stress-energy one 
obtains the following evolution equations: 



d t E a l 
d t r a 
d t A a 

d t E a/3 



d t N al3 
d t VL 



d,v a = 



d t q = 



fJe, 







Fjr p + d a q 
F a l3 A f3 + ±3 /3 £ Q/3 

(q - 2)S Q/3 - 2N <a 1 N f3> ~< + 7V 7 7 7V <a/3> 
q«* A P> + 2r <a A l3> 



d <a rP> 

~fS<a 



{d 1 -2A 1 )N !3> s + m a P 



(13) 
(14) 
(15) 



(16) 
(17) 



qN af3 + 2T, {a s N p)S - e^d^s 
{2q - 1)Q - 3P - d a Q a + 2Q a A a 
n Q/3 s Q/3 (18) 

(G-5 a p + 2v a vp){d t Q fi -2[q+l]QP) 

2v a (d t Q - 2[q + l]n) (19) 
[2{q-2) + \(2A a -r a )d a -\d a d a ]q 

2 
A 

2(q-2)tt + d t n 



\d a r a + lA a r a + |r /3 a a S a/3 - 2Y, af, W a p 



2 

2n 



2QE al3 v a vp 
v a d t v a 



(20) 



Here we are using units where c — 8ttG = 1 . Furthermore 
the quantities G_, F a p, W a p and dtQ a are given by 



G- = l-v 2 

F a /3 = qSa/3 — £q/3 

- ld a r p -\e< 5 a {d 1 -2A 1 )Np S 

d a p - d p n 



(21) 
(22) 

(23) 



dtQa = 2(q-l)Q a -J: a0 Q 13 

+ (P-0)r Q +n Q/3 (3^ + r ) 



a/3 



In addition to the evolution equations, the variables 
satisfy constraint equations as follows: 

= (C com r = e a ^{d a E p l - [r a + A a ]E fj l ) 

- N^Ej (25) 
= C G = 1 + \{28 a - 2r a - 3A a )A a ~ \N a pN^ 

+ U Na «f ''~^a0^ (26) 
= (C C ) Q = + 2r a -E 01 ^ -3A fi Z a P 

- e a ^N ps ^ S + 3Q Q (27) 



= C a 



lya/3y , la r a _ 2 \ a 



- ±(ft + 3P) (28) 

= (C. T ) a = dpN a P - (r f3 + 2A )N a ? 

+ ^(PpA^-rpA,) (29) 

= (C W ) Q = e Q/37 (V 7 - Mi) ~ N a ?r (30) 

We want a class of initial data satisfying these con- 
straints that is general enough for our purposes but sim- 
ple enough to find numerically. Recall that on an initial 
data surface, the spatial metric hij and extrinsic curva- 
ture K tJ must satisfy the constraint equations 



R + K 2 -K ij Kij = 2/z 



(31) 
(32) 



Here Di and R are respectively the derivative operator 
and scalar curvature associated with hij and fi and qi 
are the components of the stress-energy tensor given in 
equation (|3J). We use the York method^ which begins 
by defining the quantities hij and A %3 by 



hij = tp 4 hij 



K ij - \Kh l i = ^ w A lj 



(33) 
(34) 



(24) 



We choose K to be constant, hij to be the flat metric &y 
and q % to vanish. With these choices, equations l|31|) and 
(|3^|l become 

8iA ij = (35) 
= {±K 2 - i M )^ 5 - l^Aijip- 7 (36) 

Here di is the ordinary derivative with respect to Carte- 
sian coordinates and indicies are raised and lowered with 
Sij. 

We choose space to have topology T 3 with the Carte- 
sian coordinates x, y and z each going from to 2ir. We 
choose the following solution of equation (|33|) 

A 11 = 0,2 cos y + 0.3 cos z + 62 + 63 
A 22 = ai cos x — 03 cos z + b\ — 63 
A 33 — —ai cos x — a.2 cos y — 6 X — b 2 (37) 

with the off-diagonal components of A 13 vanishing. Here 
the ai and hi are constants. Note that due to the period- 
icity of the coordinates and the linearity of equation l|35|) 
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the general solution of equation l|35|l is a Fourier series. 
The solution that we choose is then essentially the sim- 
plest solution of equation (|3"5)l without symmetries. The 
quantity /i can be freely specified and we choose it to be 

/i = c + ci cos x + C2 cos y + C3 cos z (38) 

where the Cj are constants. With these choices for 
and equation (|36|) is solved numerically (in a manner 
to be described later) to yield ip and therefore and 
K ij . 

From this initial data, we must then produce the initial 
values of the scale invariant variables. From equation Q 
it follows that H = —K/3 and since H is constant it then 
follows that r a vanishes. Since the initial spatial metric 
is conformally flat, we choose the initial spatial tetrad 
vectors by 

E a l = H- V~ 2 <5a* (39) 

It then follows from equation © that N a p vanishes and 
that 

A a = -2^ -1 a V (40) 
From equation JQ) it then follows that 

while Q is given by equation JSJ and q by the vanishing 
of equation (J2SJ. 

The numerical method used is as follows: each spatial 
direction corresponds to n + 2 grid points with spacing 
dx = 2ir/n. The variables on grid points 2 to n + 1 are 
evolved using the evolution equations, while at points 1 
and n+2 periodic boundary conditions are imposed. The 
initial data is determined once equation i|3T)|) is solved. 
This is done iteratively as follows: Q Define 

Sty) ee -24> + {±K 2 - i M )V 5 - \A^A^- 7 (42) 

Then equation 1|36[) takes the form d l diip — 2tp = S(ip). 
We make an initial guess ip° for ip and solve using the 
conjugate gradient method ^3] the equation 

d l d^ k+1 - 2ip k+1 = S{ip k ) (43) 

iterating until tp k satisfies equation ll.'itill to within a set 
tolerance. 

The evolution proceeds using equations (|13I2U|) with 
the exception that the term (5 — 2q)C q is added to the 
right hand side of equation (|20|) to prevent the growth 
of constraint violating modes and the term — 0.6(Cc) is 
added to the right hand side of equation H15|) to make the 
system well posed. Spatial derivatives are evaluated 
using centered differences, and the evolution is done us- 
ing a three step iterated Crank-Nicholson methodfl^ (a 
type of predictor-corrector method) . In equation H2U|) the 
highest spatial derivative term is — ^d a d a q which gives 
this equation the form of a diffusion equation. Note that 



diffusion equations can only be evolved in one direction 
in time, in this case the negative direction which corre- 
sponds to the approach to the singularity. Stability of 
numerical evolution of diffusion equations generally re- 
quires a time step proportional to the square of the spa- 
tial step. However, the constant of proportionality de- 
pends on the coefficient of the second spatial derivative. 
To ensure stability, we define -E max to be the maximum 
value of \E a l \ (over all space and over all a and i) and 
then define dt\ = — \ {dx / E max ) and dti = —^dx. The 
time step dt is then chosen to be whichever of dt\ and 
dt2 has the smaller magnitude. 

Before presenting numerical results, it is helpful to con- 
sider what behavior to expect as the singularity is ap- 
proached (that is as t — -> —00). First denote the eigen- 
values of E Q /3 by (Si,E2,E3). Then suppose that at 
sufficiently early times the time averages of q — Sj are 
all positive. Then the time averages of the eigenvalues of 
F a j3 are all positive. Since we are evolving in the negative 
time direction, this should lead (through equation ffT^i ) 
to an exponential decrease in E a l . However, since all 
spatial derivatives appear in the equations in the form 
d a = E a l di we would expect the spatial derivatives to 
become negligible. That is, the approach to the singu- 
larity is local. Furthermore, this positivity of the time 
averages of the eigenvalues of F a p should lead (through 
equations (|14I15|) ') to exponential decrease in r a and A a . 
A similar argument applied to equations l|18|) and l|24|) 
and using equations (|26|) and l|28|l indicates that as the 
singularity is approached Q a should become negligible, 
but f2 should not, and therefore that v a should become 
negligible. Thus, as the singularity is approached, the 
system should be well described by a simplified set of 
evolution and constraint equations where spatial deriva- 
tives as well as r a , A a and v a are negligible. Note that 
the fact that spatial derivatives are becoming negligible 
does not mean that the spacetime is becoming homoge- 
neous. Rather the considerable spatial variation is be- 
coming a negligible part of the equations of motion since 
all spatial derivatives appear multiplied by E a l which is 
becoming negligible. We now write down this simplified 
system of evolution and constraint equations where all 
these terms are neglected. In this limit equation l|27|) 
implies that the matricies S a p and N a p commute and 
therefore have the same eigenvalues. The evolution equa- 
tions for S a ^ and N a @ can then be written as evolution 
equations for their eigenvalues. The non-trivial evolution 
and constraint equations then become in this limit 



9 t S, = (q - 2)Ei - 2Nf + \J2 Nk j N i + l Y ( 44 ) 

d t N, = {q + 2H i )N i (45) 
d t fl = 2(q-2)Q (46) 
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3 

= y + ^S| + 6fi-6 (47) 

fc=i 

3 

= 5>2+6f>-3 g (48) 
fe=i 

Here Sj and JVi are the eigenvalues of E Q /3 and A Q /3 
respectively and Y is given by 

fe=i \fc=i / 

and indicies are not summed over unless explicitly indi- 
cated. 

Suppose that the dynamics is in a period (called a Kas- 
ner epoch) when all the N are negligibly small. Then 
it follows from equations 147JI and H48J1 that q = 2 and 
therefore, from equations (|44|l and H46|) that and the 
Ej are constant. From equations 147|) and l|45l) it follows 
that there are two possibilities for a Kasner epoch: (i) all 
the Ej are > — 1 in which case the Ni remain negligible 
and the Kasner epoch lasts all the way to the singularity, 
(ii) one of the Ej is < — 1 in which case the correspond- 
ing Ni grows until it is large enough to bring the Kasner 
epoch to an end. We now look in more detail at possi- 
bility (ii). Let Ei be the E, that is < -1. Then Ni is 
the N that is growing. We are therefore led to examine 
equations (|44I48|) neglecting N2 and N3 but not N\. In 
this regime equations l|44|l and (|45|l become 



5 t Ex = -S 2 (E!+4) (50) 

S t E 2 = -S 2 (E 2 -2) (51) 

5 t E 3 = -S 2 (E 3 -2) (52) 

d t n = -2S 2 n (53) 



where S = Ni/y/6. Define Z = 4 + Ei. Then equation 
(|50|l becomes dtZ = —S 2 Z while from equations l|51|) and 
(|52|l it follows that there are constants c 2 and C3 with 
c 2 + C3 = —1 such that E 2 = 2 + c^Z and E3 = 2 + C3Z. 
Similarly, it follows from equation l|53|) that there is a 
constant C4 such that Q — C4Z 2 . Finally, it then follows 
from equation i|47|) that Z satisfies the equation of motion 

d t Z =[\{A-a 2 )Z 2 -AZ + &}Z (54) 

where a 2 = 1 — [12c4 + (c 2 — C3) 2 ]. Note that the quan- 
tity in square brackets vanishes at Z + and Z_ where 
Z± = 6/(2 q= a). Therefore the dynamics is a "bounce" 
from a Kasner epoch corresponding to Z_ to one cor- 
responding to Z+. Use a minus subscript to denote 
a quantity before the bounce and a plus subscript to 
denote a quantity after the bounce. We then have 
f2_|_/Q_ = (Z + /Z_) 2 . However, from the definition of 
Z it follows that Z_ = 4 + Ei^. while from the definition 
of Z± it follows that (1/Z+) + {1/Z-) = 2/3. We then 



find the following "bounce rule" relating a quantity after 
the bounce to quantities before the bounce. 

fi +=Kdb) 2 (55) 

Note that from the bounce rule it follows that f2 increases 
at each bounce. Furthermore, it follows from equation 
(|47|l (and from the fact that E Q ^ is trace-free) that the 
minimum possible value for a during a Kasner epoch 
is — 2yl — fl and therefore that no further bounces can 
happen once fl > 3/4 (though bounces may cease at 
lower values of fl). Thus we expect that at each spatial 
point there is a last bounce followed by a Kasner epoch 
that lasts all the way to the singularity. In other words, 
we expect the approach to the singularity to be of the 
Andersson and Rendall class. 



III. RESULTS 

All runs were done in double precision on a SunBlade 
2000 with n = 50 (except for a convergence test which 
also used n — 25). The equations were evolved from t = 
to t = —90. For the initial data, the trace of the extrinsic 
curvature was —1 corresponding to an initial value of 1/3 
for H. The constants a^, bi and Ci characterizing the 
initial data were 

a,i = (0.2,0.1,0.04) (56) 
b % = (1.7,0.1,0) (57) 
Ci = (0.005,0.005,0.005) (58) 

and the constant Co was 0.02. 

We would like to know whether E a l , r Q , A a and v a be- 
come negligible as the singularity is approached. In figure 
1 are plotted the maximum values (over all space, a and 
i) oflnl-Ec/l, ln|r Q |, ln|A"| andln|w a |. Note that after a 
certain amount of evolution, all these quantities steadily 
decrease. This indicates that after a certain amount of 
time spatial derivatives become negligible in the eqau- 
tions of motion and that the approximation considered 
at the end of the previous section becomes valid. 

It then follows that the interesting part of the dynamics 
is the development of the variables at a single point as a 
function of time. We now present data of that form. The 
behavior at the spatial point chosen is typical. 

The results of a convergence test are plotted in figure 
2. Here what is plotted is 4C q for the n = 50 run (solid 
line) and C q for the n = 25 run (dotted line) . Both quan- 
tities are plotted vs —t. Note that the two curves roughly 
agree in magnitude but become out of sync in time. This 
indicates second order convergence but with the system 
having sensitive dependence on initial conditions. Simi- 
lar results were obtained for the other constraints. 

Figures 3 and 4 show respectively the diagonal com- 
ponents of E a( g and N a p in the asymptotic frame, i.e. 
the frame of the eigenvectors that E"^ has at the end of 
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FIG. 1: maximum values of In \E a z \ (solid line), In |r a | (dot- 
ted line) ln|A a | (dot-dashed line) and In |u a (dashed line) vs 
-t 



FIG. 3: components of E a) 3 in the asymptotic frame vs —t 
Ei (solid line), E2 (dotted line) and E3 (dot-dashed line) 
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FIG. 2: 4Cq vs —t for the n — 50 run (solid line) and C q vs 
—t for the n = 25 run (dotted line) 
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FIG. 4: components of N a p in the asymptotic frame vs —t 
Ni (solid line), N2 (dotted line) and N3 (dot-dashed line) 



IV. CONCLUSIONS 



the evolution. For that part of the evolution where the 
approximation made at the end of the previous section 
is valid, these diagonal components are the eigenvalues 
of E Q () and N a p respectively. Note that the dynamics of 
the eigenvalues of S Q /3 consists of epochs where they are 
apporoximately constant (Kasner epochs) punctuated by 
short bounces. Furthermore the components of N a p are 
negligible except at the bounces. Also note that there 
is a last bounce and that this coincides with the most 
negative eigenvalue of becoming greater than — 1. 

In figure 5 is plotted vs —t. Note that the behavior 
of f2 is a sequence of constant values that are punctuated 
by short bounces and that the bounces in f2 occur at the 
same times as the bounces in E aj g. The sequence of the 
values of is 0.05956, 0.1607, 0.4139, 0.6231, while the 
corresponding values of the most negative eigenvalues of 
S"^ are -1.583, -1.565, -1.277, -0.6596. These values obey 
the "bounce rule" of equation l|55|l . 



These simulations support the expected picture for the 
approach to the generic singularity in a spacetime where 
the matter is a stiff fluid. As the singularity is approached 
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FIG. 5: Q, vs -t 
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the terms in the equations of motion involving spatial 
derivatives become negligible. The dynamics at each spa- 
tial point consists of a sequence of Kasner epochs punc- 
tuated by short bounces. The sequence of values of 
obeys the expected bounce rule. There is a last bounce, 
after which the dynamics is described by a single Kas- 
ner epoch all the way to the singularity, thus yielding a 
spacetime in the class of reference 0. 

What remains to be done is to treat the approach to 
the singularity for non-stiff fluids. Here the behavior 
that is expected is quite different. The BKL conjecture 
is that the matter will become negligible and the dy- 
namics of the gravitational variables as the singularity is 
approached will be that of vacuum spacetimes. The for- 
malism of 7] includes a class of non-stiff fluids, and the 
resulting equations are similar to those of the stiff fluid 



case. Nonetheless, the numerical methods of this paper 
are not adequate to treat the case of non-stiff fluids. That 
is because in non-stiff fluids shock waves form, while the 
numerical methods of the present paper are appropriate 
for smooth solutions. A shock capturing method would 
be appropriate to treat the non-stiff fluid case. 



V. ACKNOWLEDGEMENTS 

This work was partially supported by a grant from 
the National Science and Engineering Research Council 
of Canada and by NSF grant PHY-0456655 to Oakland 
University. 



[1] for a review see A. Rendall, "Theorems on existence and 
global dynamics for the Einstein equations" Living Re- 
views in Relativity (2002-6) 

[2] for a review see B. Berger, "Numerical Approaches to 
spacetime singularities" Living Reviews in Relativity 
(2002-1) 

[3] L. Andersson and A. Rendall, Commun. Math. Phys. 

218, 479 (2001) 
[4] D. Garfinkle, Phys. Rev. D65, 044029 (2002) 
[5] D. Garfinkle, Phys. Rev. Lett. 93, 161101 (2004) 
[6] V. Belinskii, I. Khalatnikov and E. Lifschitz, Sov. Phys. 

Usp. 13, 745 (1971) 



[7] C. Uggla, H. van Elst, J. Wainwright and G.F.R. Ellis, 

Phys. Rev. D68, 103502 (2003) 
[8] J.W. York, Phys. Rev. Lett. 26, 1656 (1971) 
[9] J. Isenberg, Class. Quant. Grav. 12, 2249 (1995) 
[10] W. Press, S. Teukolsky, W. Vetterling and B. Flannery, 
Numerical Recipes in FORTRAN, second edition (Cam- 
bridge University Press, Cambridge, 1992) 
[11] D. Garfinkle and C. Gundlach, gr-qc/0501031 
[12] M. Choptuik, in Deterministic Chaos in General Relativ- 
ity, edited by D. Hobill, A. Burd and A. Coley (Plenum, 
New York, 1994), pp. 155-175 



